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A METHOD FOR CALCULATING LAMINAR VISCOUS, 


COMPRESSIBLE FLOWS WITH SMALL 
PRESSURE GRADIENTS 

E. Dale Martin 
Ames Research Center 


SUMMARY 


A method is presented for calculating laminar viscous, compressible flows. Pressure gradients 
are assumed to be negligible in energy conservation and in developing property relations, but not in 
momentum conservation. 

The equations are developed from the Navier-Stokes equations in a noninertial reference 
frame. The method is an extension of an approximate method developed previously for calculating 
combined forced and contained natural convection in a rotating tank. The restriction to small 
density and temperature variations in the previous method is removed. 

For illustration, a two-dimensional simulation of the rotating tank problem is formulated, and 
the numerical computation procedure to be used is indicated. 


INTRODUCTION 


This report presents a method for calculating laminar viscous, compressible flows in which 
pressure gradients are small enough to be neglected in the governing energy equation and in certain 
fluid-property relations. Pressure gradients as sources of momentum change need not be negligible. 

The present method was developed as an extension of a previous approximate method for 
calculating combined forced and contained natural convection in a tank with time-dependent 
rotation. For use in studying the convection processes in the cryogenic oxygen tanks in a rotating 
spacecraft (see ref. 1 ), the previous method was developed in reference 2 from the Navier-Stokes 
equations. The corresponding computational method is given in reference 3 and was used with 
evaluations of thermodynamic properties from reference 4 and special methods for analyzing 
thermodynamic states and simulating an internal heater from reference 5 to obtain results described 
in reference 6 for a square-tank configuration. Subsequent corresponding results for a circular-tank 
configuration are given in reference 7. 

The approximate theory of reference 2 represents a small-density-variation approximation of 
the Navier-Stokes equations that is an extension of the Boussinesq approximation to account for all 
effects of time-dependent rotation. It includes the influence of the combined effects both of 
effective buoyancy body forces due to temperature and density stratifications (natural convection) 



and of arbitrary time-dependent rotation of the container (forced convection). That method is 
restricted to small density and temperature variations and was considered adequate for its intended 
application (ref. 8): a study of convection in supercritical cryogenic oxygen in a gravitationless 
environment (with only very small body forces due to rotation). 

The present method removes the restrictions of small temperature and density variations. It 
accounts for important effects of compressibility (except for sound waves) and variable properties 
in viscous flow. It can be regarded as a further generalization of the Boussinesq approximation, but 
is not restricted to either natural convection or convection in a closed container. Since the primary 
assumption in the method (besides laminar flow) is that pressure gradients be sufficiently small, the 
method should apply to any laminar flow regions in which that assumption is true (flow regions in 
which density changes are due mainly to changes in temperature distribution and to changes in 
average pressure level). This includes convection in a rotating tank with rather large density and 
temperature differences and, for example, some regions of separated flow over bodies (see ref. 9 for 
a comprehensive survey). 

As in reference 2, all relevant effects of variable rotation are included in the full 
three-dimensional equations. The final equations and conditions are in a form similar to those of 
reference 2, but with significant differences. A stream function is not conveniently used in this 
method, but the Poisson equation for the stream function in reference 2 is just replaced by 
corresponding Poisson equations for the velocity components for use in one step of the 
computation procedure. The rotating-tank problem is formulated in a form similar to that in 
reference 2, and it is expected that the numerical methods of reference 3, as well as the special 
procedures of references 4 and 5, could be easily adapted to this formulation. The computation 
procedure is briefly sketched. 


NAVIER-STOKES EQUATIONS IN A TIME-VARYING 
ROTATING REFERENCE FRAME 


The basic differential equations governing compressible flow of a Newtonian fluid in a 
noninertial system are the Navier-Stokes equations in the form (see ref. 2): 
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V= -\.DR 

~ 9 Dt 


a) 


dv , , 

p ^7 = -Vp+p^ +V* [X(v • V) n +v -|p[VF+(VF) ? ] > (2) 
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where: D/Dt is the substantial-derivative operator, 


_D_ 

Dt 
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dt 


+ V • V 
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the apparent body force per unit mass in the noninertial system is 

g = -[a + Q X R + 212 X V + Q X (g X £)] (4) 

the viscous dissipation of energy is 

p<F = X(V • F) / ; VK + m[VF+(VK)^]:VF (5) 

and 

F fluid velocity in the noninertial system 

p local mass density 

p local thermodynamic pressure 

T temperature 

p shear viscosity coefficient 

X k — (2/3)p is the second viscosity coefficient 

k bulk viscosity coefficient 

1^ unit tensor 

(VF)j transpose of the dyadic vF 

h specific enthalpy 

k thermal conductivity in the Fourier heat conduction law for heat flux 

£ (possibly time-dependent) linear acceleration of the noninertial system (in excess of 

SX X R) relative to an inertial system 

12 time-dependent angular velocity of the noninertial system relative to an inertial system 

R perpendicular vector from the axis of rotation to a point of interest in the flow system 

(see ref. 2) 

£1 angular acceleration of the noninertial system, dSl/dt 

• • 

In equation (4), Q X R is equal to a linear acceleration due to the angular acceleration Q, 2H X V is 

the Coriolis acceleration, and Q X (Q X R) is equal to the centripetal acceleration due to rotation 
of the system. The direction of the axis of rotation is assumed to be constant— 

Q/£l = £ 3 (6a) 
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where £2 is the time-dependent magnitude of £2 and <? 3 is a constant unit vector in a direction 
parallel to £2. Note that by the definition of R 

£2 • R = 0 (6b) 

/-V/ />-< 

To make the system of equations complete, one must specify auxiliary expressions for the 
transport properties p, X, and k, and equations of state, 

P = p(p,T) (7a) 

h = h{p,T) (7b) 


ASSUMPTIONS AND RESULTING VECTOR EQUATIONS 


In flows in which pressure gradients are sufficiently small so that density changes are due 
mainly to changes in temperature distribution and to changes in the average pressure level at a given 
time, the terms V • Vp and p4> can be omitted from the energy equation (3). Those terms are of 
minor importance in that equation whose primary function is to determine the changes in the 
temperature distribution. Refer also to Whitham (ref. 10, pp. 124-127), who shows that those 
terms are negligible if 


T-T j 
T i 


»M l 2 


where T x is a typical temperature in the flow andA^ is a corresponding typical Mach number. The 
method of reference 2 also required {T— T 1 )/T 1 to be small; the present method does not. The 
above conditions allow the pressure p to be replaced by its instantaneous spatial average p(t) for use 
in the energy equation: 


h(p,T)**h[m, T ] 

QE ss =p'(t) 
Dt dt P W 


Similarly, whereas in reference 2 use was made of 


P(P, T) - p Q + (T- T 0 ) 


(If) 


+ 0 [(AT) 2 ,Ap] 


PJO 


( 8 ) 

(9) 


we write now the more accurate representation based on the above conditions: 

p(p,T) = p(p 1 ,r) + (p-p 1 )(^^+0[(p-p 1 ) 2 ] ^p{p x ,T)+[p{t)- Pl ](^i T (10a) 


In the above, AT =T — T 0 and A p-p — p Q , where p Q , p Q , and T 0 in reference 2 represented the 
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spatial average thermodynamic state; and p x is a constant typical pressure in the flow. Note that 
equation (10a) includes an assumption of small pressure gradients because p is replaced by pit). 

Similarly, the transport properties p and k, as well as the thermodynamic state properties 
including (dh/dp)p and 


C P 



can be approximated in the form (cf. (10a)) 


X(p,T) = F l (T) + 


p(0~pA 

k P2 ~Pi / 


F 2 (T) 


( 11 ) 


(10b) 


where p x and p 2 are constants and F x (T) and F 2 (T) are different functions to be specified for each 
property X to be so represented. The validity of this approximation was found from experience in 
the progress of studies reported in references 5 and 6. With h in (7b) represented approximately by 
equation (8) and with the definition (11) and 


c t = 



( 12 ) 


one obtains from (8): 


Dh = 
Dt 


DT. 

C PDt +Ct 


(13) 


In addition to the general property representation (10b), it is convenient to write specifically (for 
(10a)) 


p = p 0 [F x {T)+f{t)F{T)] 


(14) 


in place of (7a), where p 0 is a constant, F x (7) and F{T) now denote definite functions to be 
specified, and fit) is the bracketed factor in (10b). 

With the above assumptions the energy equation (3) yields 


= — V • ikVT) - — + — 
Dt pCp ~ ~ Cp pCp 


Pit) 


(15) 


It is then convenient to define 


a — — 


LM 

P Dt 


(16) 
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and to use equation (14) to obtain 


« = - J-{f\OF(T) + [ f{t)F\T ) + F\ (D] 


(17a) 


where 

( )' generally denotes the derivative of a function 

p a /p is obtained from (14) 

DT/Dt is specified by (1 5) 

Therefore 


a= ~^r{f'(t)F(T) + [f(t)F'(T)+F[m] 


— v • (kVT) - — + -L 

pCp ~ ~ Cp pCp 



in which p is given by ( 1 4). 


(17b) 


The resulting vector equations are now, from (1) and (16), conservation of mass: 


(18) 

where a is given by (17b); conservation of momentum given by equation (2); and from equations 
(15) and (18) a “divergence form” (convenient for numerical computation) of the energy 
conservation equation: 


^ + V • {VT) = aT + -^~ V (kVT )~— + — p'(t) 
9 1 ~ ~ pc p ~ c p pCp “ 


(19) 


FORMULATION IN TERMS OF VORTICITY 


To use a solution procedure similar to that described in reference 3, it is convenient to 
formulate the problem in terms of vorticity. In addition to the appropriateness of vorticity as a 
primary variable (see Lighthill, ref. 11, pp. 57-60), it is convenient to replace equation (18) by a 
corresponding vector Poisson equation. The formulation also eliminates from the problem both the 
second viscosity coefficient X and the pressure, except for the average pressure p(t), which must be 
determined from an appropriate equation of state as in reference 5 in terms of a spatial average 
temperature T and corresponding average density p. 

With the definition of vorticity 

co = V X V (20) 
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and the identity 


VXu-VX(\7XF) = V(VF)-V 2 ^ (21) 

/-s^ /-N./ /-s^ /"V, /-s^ 

along with (18), one obtains a vector Poisson equation for V in terms of oj and a: 

V 2 V=va~ vX gj (22) 

/-s»/ 

which would be solved for the velocity components after co is determined at each time step. 

The momentum equation (2) is put into the form of a vorticity equation as follows: Use the 
identity 


V- -VX (VX V) 


(23) 


to write equation (2) as 

= V(-p + Xa) + + V • |ju[ VV + (ViO/]| 

Using the identities (for any scalar 0 and vector A ): 

V X (4>A) = 0V X A + V0 X A 
V X V0 = 0 


dV 

dt 


- 


FX co 


(24) 


(25a) 

(25b) 


take the curl of equation (24) to obtain 


0GJ 

Jt 


-VX (FX co) 


DV c \ 

Vp X ^ = pv x ^ + VP x ^ + V x V • <ju[ VF + (VF)f] i 


(26) 


For use in equation (26) 


1 [f(.t)F\T) + F\{T) ] 

- vp = -= vr 

[f(t)F(T) + F l (T)\ ~ 


(27) 


The following developments are also useful for reducing equation (26). From equation (4) along 
with 


n X (ft X R) = -FL 2 R = - v(^n 2 i? 2 ] (28) 

£ = V0 (29) 

(where <p = R* ' £(t) is a scalar potential for the time-dependent acceleration £ and is the 
position vector from an arbitrary fixed point on the axis of rotation in the noninertial system; see 
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ref. 2), one can write 


£ = ~ 



±n 2 R 2 ) + &XR + 2£2X V 

Z f ~ ~ ~ ~ 


(30) 


Then 


V X g = -[ V X (Cl X R) + 2V X (Cl X 7)1 
For arbitrary vectors ^4 and Bj the identity 

V X 04 XB) = B -VA-A ' VB + AV - B - BV- A 

/^/ /■*+/ y^yy^y r**y 


= V (BA)-V’ (AB) 


(31) 


(32a) 

(32b) 


yields 

VX(ilXR) = Clv R = 2CI (33a) 

and 

V X (ft X F) = -n • vF + fia (33b) 

y^ y^y y^y y^yf^y /*■»✓ 

= -Cl (dV/dZ) + Qa (33c) 

where Z is a coordinate in the direction of the unit vector e 3 defined by equation (5). Therefore, 
results for use in equation (26) are 


g = ~(a + Q X R + 2Q X V - Cl 2 R) 
dV 

V X g = 2(Cl -Qa-Q) 


(34a) 

(34b) 


Further, if the identity (32b) is used in the left side of (26) for VX(KXw), the above 
developments yield 



where 

^-Vp is given by (27) 

H 

g is given by (34a) 
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DV/Dt =dV/d t+V-VV 

p is given by ( 1 4) 

The last term in equation (35) can be reduced further for possible convenience in the solution 
procedure. The identities: 

V * [/uVF] = pV 2 V + (VP) ‘VV (36a) 

X ' [m(VK)^] = pV(V * V) + (Vp) • (Vf> (36b) 

can be used along with equation (22) to obtain 

V ' |m[VF + (VF)^ | = M[2Va- VX o>] + (Vju) * [VK + (ylO ? ] (37) 

When the curl is taken of equation (37) for use in (35), one can use the identities (25) along with 
the additional identities: 

D = ° ( 38a) 

VX (VX C 0 )=V(V- CO) -v 2 CO ^-v 2 CO (38b) 

and with equation (22) to obtain for the last term in (35): 

p V x + = p S + <S“> x <2a+v 2 F)j 

+ j2 x {(^>Mvr+(vF),]} (39) 

The final results for this formulation are now: 

1. Equation (14) for p, where f(t) is the bracketed factor in (10b) 

2. Equation (17b) for the auxiliary variable a 

3. The vector Poisson equation (22) for V 

4. Equation (35) for co, where the last term may be replaced' by equation (39), (l/p)Vp is 
given by (27), and g is given by (34a) 

5. Equation (1 9) for T 

6. Equation (1 2) for c t 

These equations must be supplemented by appropriate relations in the form of (10b) for Cp , 
( dh/dp p , and k , and appropriate relations for p = p(t) as a function of p and T and for F{T) and 
Fi(T) in equation (14). 
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Figure 1 . — Square tank configuration. 
The vorticity vector is 


EQUATIONS AND CONDITIONS FOR TWO- 
DIMENSIONAL CONVECTION IN A SQUARE TANK 


For a two-dimensional simulation of the convection 
in a tank it is convenient to use a square-tank configura- 
tion with the equations and conditions in Cartesian 
coordinates. For convection in a plane normal to the axis 
of rotation, the assumed configuration is as shown in 
figure 1. With R c being the distance from the center of 
rotation to the tank center and R the radius vector from 
the center of rotation, the components of R in terms of 
the tank-fixed coordinates (x, y) are 

Ri = x - (40a) 

R 2 =y-^l + R c (40b) 


cu = oj 


where 


, , _ 3v _ 9w 
3x 3 y 


(41a) 


(41b) 


and where u and v are the x and y components of V. The x and y components of the apparent body 
force per unit mass j* are, from (34a), 


gi = R : £2 2 + 2£2v + R 2 £l ~a x (42a) 

g 2 - R 2 £l 2 ~ 2£2w - a 2 (42b) 

where a x and a 2 are the components of £. 

For two dimensions we now have, from (14) and (17b): 


p = p 0 [F iiT) + f {t)F(T)\ 


(43) 


and 


a = 


- p o h> 


f'(t)F(T ) + [f{t)F'{T) + FKD] 




2i 

c p 



( 44 ) 
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Equation (22) yields two scalar Poisson equations: 


3 2 w 

+ a 2 ^ . 

_ 3a . 

_ 9a) 

3x 2 

3y 2 

3a: 

dy 

3 2 v 

+ 3 2 v . 

. 3a . 

f 3cp 

3 a: 2 

3y 2 

dy 

dx 


Denote by A in equation (27) the factor 


f(t)F\T) + F[{T) 
f(t)F(T) + F 1 (T) 


(45a) 

(45b) 


(46) 


and by W the factor multiplying 1/p in the viscous diffusion term (the last term) in equation (35). 
Then one obtains from (39): 


W = n(¥“ +^\ + *±($* + ¥1 + \ 

V 3 a 2 dy 2 ) dx V d y dx 2 dy 2 1 

_ dj^ /9a , 3 2 m , 3 2 m \ . 3 |~ 3p , / 3v_ + + /? \~| 

a y\dx dx 2 dy 2 ) ax|_ 3 * [dx dy) dyVdy)\ 

_ 9 TSp A du\, dp (dv , du \~| 
dy l_3x \ dx) dy \ dx dy )\ 


(47a) 


or more directly from the left side of (39), 


W = 


3 2 



+ 92 [2 

p(*L- 

3m \ 

-iir 


, 

3a 2 

L M \ 3a 

3T/J 

3 a 3y L 


3a )_ 

3y 2 L 

\ 3a 

dy )_ 


(47b) 


In terms of A and W the vorticity equation (35) becomes in two dimensions: 


Jf( 


g 2 


Dy\_dT 

Dt 


)-&(«■ 


Du \ 
Dt) 


-2(S2a + f2) + -|-lV (48) 

r 


where 


Du _du_ + du +y du 
Dt dt dx d y 


+V |L 

Dt dt 3a dy 


(49a) 

(49b) 


The two-dimensional energy equation for T with the convective terms in the divergence form is, 
from (19), 


|f + A (uT) + |- 
3r 3 a dy 


(vT) = aT+-!~ 


_ _3_ 

/*»r 

Wi 


_3a 

r M , 

I dy 1 

1 sWJ 


+ 


c p P c p 


p\t ) 


(50) 
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In equations (44) and (50), c t is given by equation (12), and 


m = 


Pit)-p i 
Pi ~Pi 


( 51 ) 


(See eq. (10b)), where p(t) can be computed by a method described in reference 5, from a 
knowledge of the mean fluid density p and the temperature distribution. 


Initially one can specify the conditions, 

/ = 0 : u = v - go = 0) 

' all x,y 

= £ = °\ < 52 > 

if the calculation is started at an initial condition of “rigid body rotation” before the rotation rate is 
changed. Also an initial temperature distribution can be specified. One must specify 12(7) at all 
times, including t = 0. 


At all boundaries, the conditions of no slip in the tangential direction and no normal mass flow 
are represented by the conditions 


u = v = 0 on all boundaries (53) 

for application to the Poisson equations (45). Boundary conditions on co needed for application to 
the vorticity equation (48) (with either (47a) or (47b) for the viscous diffusion) are supplied as 
follows. If subscript b denotes a value on the boundary, subscript Ax denotes a value at a small 
distance |Ax I inside the boundaries x = 0 and x = / (see fig. 1), and subscript A y denotes a value at a 
small distance |A y I inside the boundaries y = 0 and y = /, then from the definition of vorticity (41b) 
with the conditions (53), one has 


a 1 *' 0 = 

at y = 0,1: = 


which can be represented for numerical approximation by: 

0)Ax 


x = 0 c op - 


| Ax I 


X = l CO/y 


_ ~0)Ax 
|Ax | 


(54a) 

(54b) 


(55a) 

(55b) 
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(55c) 


7 = 0 


« 'b 


( u )Ay 

~WT 


y = i 


( w ) Ay 

^ = WT 


(55d) 


For application to equation (50), one may specify dT/dx on x = 0 and x = l and dT/dy on y = 0 and 
y = /. In addition, a heat source at some interior point can be specified to represent the power input 
from a heater (see refs. 1 and 5). 


SOLUTION PROCEDURE 


For convenience, the above equations, initial conditions, and boundary conditions could be 
put into dimensionless form as was done in reference 2 for the approximate method. It is 
anticipated that the numerical solution procedure of reference 3 for the approximate equations 
from reference 2 could be adapted to the present formulation. The dimensionless differential 
equations corresponding to the equations in the above section would be converted to difference 
equations in much the same manner. One new aspect is the appearance of Du/Dt and Dv/Dt in 
equation (48) and of p(t) and p f ( t ) in several equations. The time derivatives du/dt, dv/dt, and 
p f ( t ) in those terms are expected to be small and can presumably be backward-differenced without 
significantly affecting the numerical stability, because the vorticity and energy equations are being 
solved for co and T at each time step with u and v determined by Poisson equations and p from a 
state relation. Another new aspect is the replacement of the Poisson equation for the stream 
function in references 2 and 3 by two Poisson equations for u and v that depend on space 
derivatives of co at each time step. The latter two equations must satisfy both the no-slip and the 
no-mass-flux boundary conditions, which means that the vorticity at the advanced time step is not 
required to satisfy any independent conditions as in the earlier procedure (cf. also the discussion of 
boundary conditions in ref. 7). 

The suggested computational procedure for the corresponding difference equations is: 

1. Set the variables initially according to equations (52) and the prescribed temperature 
distribution. Also set p (0) at the desired initial pressure and p (0) = 0. 

2. Calculate the “predicted” vorticity and temperature distributions by the predictor in 
MacCormack’s explicit method (ref. 12; cf. ref. 3) from the difference equations corre- 
sponding to equations (48) and (50), and using conditions (55). 

3. Calculate the corresponding p by a method described in reference 5 using the temperature 
distribution and the known mean fluid density; then calculate p and a according to 
equations (43) and (44). 

4. Calculate the corresponding “predicted” u and v from equations (45) using Buneman’s 
Poisson solver (ref. 13; cf. ref. 3) with conditions (53). 

5. Calculate the corrected vorticity and temperature distributions using the corrector in 
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MacCormack’s method (cf. ref. 3) for the difference equations from (48) and (50) (and 
using conditions (55) with predicted velocities). 

6. Calculate the corrected p, p, and a. as in step 3. 

7. Calculate the corrected u and v as in step 4. 

8. Advance the time index and repeat steps 2 through 8. 

The procedure is continued for the number of time steps desired. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, June 21, 1972 
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